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ABSTRACT 

Hierarchical two-planet systems, in which the inner body's semi-major axis is between 0.1 and 0.5 AU, usually present high eccen- 
tricity values, at least for one of the orbits. As a result of the formation process, one may expect that planetary systems with high 
eccentricities also have high mutual inclinations. However, here we show that tidal effects combined with gravitational interactions 
damp the initial mutual inclination to modest values in timescales that are shorter than the age of the system. This effect is not a direct 
consequence of tides on the orbits, but it results from a secular forcing of the inner planet's flattening. We then conclude that these 
hierarchical planetary systems are unlikely to present very high mutual inclinations, at least as long as the orbits remain outside the 
Lidov-Kozai libration areas. The present study can also be extended to systems of binary stars and to planet-satellite systems. 
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1. Introduction 

Today, nearly 100 multi -planet systems have been reported, 
of which roughly 1/3 possess at least one "moderate close-in 
planet", that is, a planet with a semi-major axis between 0. 1 and 
0.5 AU 1 . Planets in this range are supposed to undergo signif- 
icant tidal interactions, resulting in slowly modified spins and 
orbits. However, for the typically assumed dissipation rates for 
gaseous planets, the spin of moderate close-in planets reaches 
an equilibrium state in only a few million years, while the or- 
bital evolution can last for the entire age of the system (Gyr 
timescale). 

Among two-planet systems, there is a special class whose 
semi-major ratio a\jai is lower than 0.1, the so-called "hierar- 
chical systems". This class counts at least 20 members (1/5 of 
all multi-planet systems), and usually at least one of the planets' 
orbits is highly eccentric 1 . During the formation process, the or- 
bital eccentricities can increase through gravitational scattering 
(e.g. Nagasawa et al. 2008), which is simultaneously responsible 
for an increase of the orbits' mutual inclination (e.g. Chatterjee 
et al. 2008). Evidently, smaller mass planets, that are as yet un- 
detected, can exist in these systems, but for a semi-major axis 
larger than 0. 1 AU the orbits usually present high values in the 
eccentricities, which reduces the stability areas for additional 
companions. 

Hierarchical systems are particularly interesting from a dy- 
namical point of view, because they can be stable with very ec- 
centric and inclined orbits, which is an uncommon behavior. In 
particular, they are very interesting when the inner planet is suf- 
ficiently close to the star to undergo tidal interactions, since the 
final outcome of the evolution can be in a configuration that is to- 
tally different from the initial one. Because the semi-major ratio 
is small, low-order mean motion resonances cannot occur, which 
allows us to perform analytical approximations such as averag- 
ing the orbits over the mean anomalies. In addition, tidal effects 
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usually act over very long timescales and therefore approximate 
theories also allow one to accelerate the numerical simulations 
and to explore the parameter space much more rapidly. 

Secular perturbation theories based on series expansions 
have been used for hierarchical systems. For low eccentricity 
values, the expansion of the perturbation in eccentricity series 
is very efficient (e.g. Wu & Goldreich 2002), but this method is 
not appropriate for orbits that become very eccentric. An expan- 
sion in the ratio of the semi-major axis 01/02 is then preferred, 
because exact expressions can be computed for the secular sys- 
tem (e.g. Laskar & Boue 2010). The development to the second 
order in ai/a 2 , called the quadrupole approximation, was used 
by Lidov (1961, 1962) and Kozai (1962) for the restricted in- 
ner problem (the outer orbit is unperturbed). In this case, the 
conservation of the normal component of the angular momen- 
tum enables the inner orbit to periodically exchange its eccen- 
tricity with inclination (the so-called Lidov-Kozai mechanism). 
For planar problems, the series expansions in 01/02 should be 
conducted to the octupole order (e.g. Marchal 1990; Ford et al. 
2000; Laskar & Boue 2010), because the quadrupole approxi- 
mation fails to reproduce the eccentricity oscillations (e.g. Lee 
& Peale 2003). Since we do not have any restrictions for the ec- 
centricities or for the mutual inclination, we need to expand the 
gravitational potential in ai /a2 to the octupole order. 

The ultimate stage for tidal evolution is the spin synchro- 
nization and orbital circularization (e.g. Correia 2009). Indeed, 
the observed mean eccentricity for planets and binary stars with 
a 1 < 0.1 AU is close to zero within the observational limitations 
(e.g. Pont et al. 201 1). Although tidal effects modify the spin on a 
much shorter timescale than they modify the orbit, synchronous 
rotation can only occur when the eccentricity is very close to 
zero: the rotation rate tends to be locked with the orbital speed 
at the periapsis, because tidal effects are stronger when the two 
bodies are closer to each other. In addition, if there is a compan- 
ion body, the eccentricity oscillates (e.g. Mardling 2007; Laskar 
et al. 2012), and the rotation rate of the planet shows variations 
that follow the eccentricity. This is exactly what is observed for 
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Mercury, whose average orbital eccentricity is around 0.2, and 
its rotation is captured in a 3/2 spin-orbit resonance (Correia & 
Laskar 2004, 2012). Therefore, spin and orbital evolution cannot 
be dissociated, and some unexpected behavior can be observed, 
such as a secular increase of the eccentricity (e.g. Correia et al. 
2011,2012). 

In this paper we intend to intensify the study of hierarchical 
two-planet systems, in which the inner orbit undergoes tidal dis- 
sipation. We present here another counterintuitive behavior, the 
inclination damping, which is also a consequence of the above- 
mentioned eccentricity pumping when the two orbits are not 
coplanar. 

2. Model 

We considered a system consisting of a central star of mass mo, 
an inner planet of mass m\, and an outer companion of mass m 2 . 
We used Jacobi canonical coordinates, with ri being the posi- 
tion of m\ relative to mo, and r 2 the position of m 2 relative to 
the center of mass of m\ and mo- We assumed that the system 
is hierarchical, thus |ri| <sc |i*2l- For simplicity, we express all 
the angles in the invariable plane of the system, i.e., the plane 
perpendicular to the total angular momentum 



H = L + Gi +G 2 



(1) 



where L is the rotational angular momentum of the inner planet, 
and G, the orbital angular momentum of each body. 

The inner planet is considered an oblate ellipsoid with grav- 
ity field coefficients given by 7 2 , rotating about the axis of max- 
imal inertia (gyroscopic approximation), with rotation rate a>, 
such that (e.g. Lambeck 1988) 



J 2 = k 2 



CO 



30m ! 



(2) 



@ is the gravitational constant, R is the radius of the planet, and 
k2 is the second Love number for potential (pertaining to a per- 
fectly fluid body). We furthermore assumed that the obliquity of 
the planet to the orbital plane is zero (s = 0°), that is, L and Gi 
are aligned. Therefore, the angle between the two orbital planes, 
i.e., the mutual inclination /, satisfies the relation 



H 2 = (L + Gi) 2 + G\ + 2(L + Gi)G 2 cos / . 



2.1. Conservative motion 



(3) 



Because we are interested in the secular behavior of the system, 
we averaged the equations of motion over the mean anomalies 
of both orbits. In the invariable plane, the averaged potential, 
quadrupole-level for the spin (e.g. Correia & Laskar 2010a), 
octupole-level for the orbits (e.g. Ford et al. 2000; Laskar & 
Boue 2010), and with general relativity corrections (e.g. Touma 
et al. 2009) is given by 
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and (f is the angle between the directions of the periastrons: 

cos ip = - cos cox cos C02 - sin coi sin a>2 cos / . (9) 

a,- is the semi-major axis (that can also be expressed using the 
mean motion n,), e, is the eccentricity, and w,- is the argument of 
the periastron. We also have [3\ = mom\j(mo + m\), /?2 = {mo + 
m 1 )m2/(mo + nii+ m 2 ), G,- = y6,n,a 2 (l - e 2 ) 1 ^ 2 , and L = ^m\R 2 to, 
where £ is the normalized moment of inertia. 

The contributions to the orbits are easily obtained using the 
Lagrange planetary equations (e.g. Murray & Dermott 1999): 
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In addition, since the variations in e\, e 2 and / are related by the 
conservation of the total angular momentum (Eq. 3), we have 
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As we assumed that L/Gi ~ (R/ai) 2 <K 1, if we neglect first- 
order terms in L/Gu we simply have (z + j = 1,2): 
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When / = 0°, we have ip = ip\ - ip2 - 0)2 - <^i- Note also 
that the longitude of the node does not appear in the equations 
of motion (Eqs. 13-16) because we used the invariable plane as 
the reference plane (Eq.4), and thus A£2 = 180°. 

2.2. Tidal effects 

In our model, we additionally considered tidal dissipation raised 
by the central star on the inner planet. The dissipation of the me- 
chanical energy of tides in the planet's interior is responsible for 
a time delay Af between the initial perturbation and the maximal 
deformation. Because the rheology of planets is poorly known, 
(15) the exact dependence of At on the tidal frequency is unknown. 
Several models exist (for a review see Correia et al. 2003; 
Efroimsky & Williams 2009), but for simplicity we adopted a 
model with constant Af, which can be made linear (Singer 1968; 
Mignard 1979). The contributions to the equations of motion are 
given by (e.g. Correia 2009; Correia et al. 2011) 
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Table 1. Hierarchical two-planet non-resonant systems with 0.1 < a\ < 0.5 AU and a\ja-i < 0.1. 



Star 


Age 


mo 


mi 


nil 


0[ 


e\ 


«2 


ei 


¥>* 


R 


'max 


Ref. 


(name) 


(Gyr) 


(M ) 


(Mj) 


(Mj) 


(AU) 




(AU) 




(deg) 


Cftju P ) 


(deg) 




HD 190360 


7.8 


0.96 


.057 


1.50 


0.128 


0.01 


3.92 


0.36 


219. 


0.54 


~ 40* 


[1] 


HD 38529 


3.3 


1.48 


1.13 


17.6 


0.131 


0.25 


3.70 


0.36 


269. 


1.13 


~ 40* 


[2] 


HD 11964 


9.6 


1.12 


.079 


0.62 


0.229 


0.30 


3.16 


0.04 


307. 


0.60 


~ 40 


[3,4] 


HD 147018 


6.4 


0.93 


2.12 


6.56 


0.239 


0.47 


1.92 


0.13 


251. 


1.25 


~ 40 


[5] 


HD 168443 


9.8 


0.99 


7.66 


17.2 


0.293 


0.53 


2.84 


0.21 


252. 


1.51 


~ 50 


[6] 


HD 74156 


3.7 


1.24 


1.88 


8.03 


0.294 


0.64 


3.40 


0.43 


266. 


1.23 


~ 20 


[7,8,9] 


HD 163607 


8.6 


1.09 


0.77 


2.29 


0.360 


0.73 


2.42 


0.12 


186. 


1.05 


~ 30 


[10] 



Notes: All masses m, correspond to minimum values (/, = 90°), except for HD 38529, which has /, = 48°; <p* = of 2 - \ R was estimated 
using Eq.(41); The maximal inclination / max was estimated using k 2 At = 100 s, (*) for HD 190360 and HD 38529 starting with «[ = 0.2 AU and 
e, = 0.25. References: [1] Vogt et al. (2005); [2] Benedict et al. (2010); [3] Butler et al. (2006); [4] Wright et al. (2009); [5] Segransan et al. 
(2010); [6] Pilyavsky et al. (2011); [7] Naef et al. (2004); [8] Bean et al. (2008); [9] Meschiari et al. (2011); [10] Giguere et al. (2012). 
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We neglected the effect of tides over the argument of the pe- 
riastron, as well as the flattening of the central star. Their effect 
is only to add a small supplementary frequency to io\, similar 
to the contributions from the general relativity (for a complete 
model see Correia et al. 2011). Expression (28) for the inclina- 
tion is zero, because we assumed the obliquity of the planet to 
be zero (s = 0°). 

Under the effect of tides alone, the equilibrium rotation rate, 
obtained when d> = 0, is attained for (Eq. 25) 



fi(ei) 
Mei) 



1 + 6e\ + 0{e\) . 



(37) 



Usually, K' <k K, so tidal effects modify the rotation rate much 
faster than the orbit. It is thus tempting to replace the equilibrium 
rotation in expressions (26) and (27). With this simplification, 
one always obtains negative contributions for ci\ and e\ (Correia 
2009), 
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(39) 



(40) 



Thus, the semi-major axis and the eccentricity can only decrease 
until the orbit of the planet becomes circular. However, planet- 
planet interactions can produce eccentricity oscillations with a 
period shorter than, or comparable to, the damping timescale of 
the spin. In that case, expression (37) is not satisfied and multi- 
planetary systems may show non-intuitive eccentricity evolu- 
tions, such as eccentricity pumping of the inner orbit (Correia 
etal. 2012). 



g. 0.7 
a? 

<£ 0.6 
0.5 
0.4 
0.3 
0.2 
0.1 




log m (M Jup ) 

Fig. 1. Radius versus the mass of the planet. We plot all known close-in 
planets in the range 0.1 < a { < 0.5 AU, for which the radius was deter- 
mined by the transiting method. We observe that the radius decreases 
with the mass in a relatively regular way. 

Table 2. Stability analysis of the HD 74156 system for different sets 
(wi,/) of initial conditions (Fig. 2). Chaotic diffusion is present when- 
ever D > IQ- 6 (in bold). 



trajectory 




h 


2n/g 


log£> 


(color) 


(deg) 


(deg) 


(kyr) 




circulation 


red 


90 


10 


13.657 


-9.41 


blue 


90 


20 


15.065 


-9.12 


green 


90 


30 


18.554 


-7.66 


brown 


90 


38 


20.877 


-4.78 


gray 


30 


60 


16.992 


-5.90 


orange 


30 


70 


17.773 


-5.33 


purple 


30 


80 


20.289 


-1.60 


Vibration 


pink 


90 


42 


15.482 


-9.08 


mangenta 


90 


50 


9.918 


-7.20 


cyan 


90 


60 


7.703 


-8.50 


yellow 


90 


70 


6.804 


-9.99 


black 


90 


83 


6.330 


-7.60 



3. Application to exoplanets 

In the following sections we apply the model described in 
Section 2 to different configurations of hierarchical two-planet 



systems. To observe the damping effect of the mutual inclina- 
tion, the spin of the planet must be fully damped, but not its or- 
bit, that is, K' <K K (see Appendix A). In addition, the damping 
timescale of the spin should be of the same order as the period 
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90 180 270 360 90 180 270 360 




90 180 270 360 90 180 270 360 



to; (deg) coj (deg) 



Fig. 2. Possible secular trajectories for the HD 74156 system seen in the (coi,I) plane (left), and in the (tOi, e t ) plane (right). We show the trajectories 
using the quadrupolar approximation (top), corresponding to the level curves of constant energy, and the octupolar approximation (bottom). The 
dashed black curves in the (u)\,I) plane correspond to the separatrix between the circulation and libration zones of a)\. The colors are preserved 
in all pictures, each one corresponding to a given value of the total angular momentum of the system (determined by different values of the initial 
mutual inclination). All trajectories are compatible to the present knowledge for this system (Table 2). 



of the eccentricity oscillations, K ~ g (Eq. A. 18). This is valid 
for gaseous planets roughly within 0.1 < a\ < 0.5 AU, which 
we call "moderate close-in planets". In Table 1 we list all hier- 
archical systems known to date whose inner orbit satisfies the 
above condition. We focus our analysis on the HD 74156 sys- 
tem, but all the main results are easily extended to the remaining 
planetary systems. 

3. 1 . Radius of close-in exoplanets 

According to expression (A. 18), the v\ frequency (Eq. 18) is a 
key parameter for the observation of the eccentricity pumping 
of the inner orbit and consequent damping of the inclination (by 
means of g x ). The minimum masses and the semi-major axis are 
relatively well determined from the observations, so the largest 
incertitudes in vj come from the Love number ki, and particu- 
larly from the radius of the planet, which appears as a power of 5. 
Therefore, a correct estimate of the planetary radius is necessary 
to observe some effect on the inclination. 

Since the radius of the planet is correlated with its mass, one 
solution is to adopt a constant value for the density, p, and then 
compute the radius simply as R 1 = 3;«i/(4^p). However, by ap- 
plying this strategy to the two largest planets of the solar system, 
Jupiter and Saturn, we immediately see that it can give very dis- 



tinct results. The density of a planet depends on many factors, 
such as the age of the system, the initial composition of the ac- 
cretion disk, or where the planet formed in the disk. Any theoret- 
ical estimation of the radius is then subject to large incertitude, 
and only direct observations can give reliable values. 

We used an empirical expression based on the previously ob- 
served radius of close-in planets in the range 0.1 < d\ < 0.5 AU. 
We found ten planets in this range 1 whose radius were deter- 
mined by the transiting method (Fig. 1). We observe that the ra- 
dius decreases with the mass in a relatively regular way, there- 
fore we performed a linear regression of the observational data: 

R/R hlv = 0.46 log (m/M Jup ) + 1.10 . (41) 

This expression also agrees well with the solar system data, giv- 
ing 1.1 Rj up for Jupiter and 1.0/?sat for Saturn. 

3.2. Initial conditions uncertainty 

Assuming that the observational values of the minimum masses, 
semi-major axis, and eccentricities of the planetary systems 
listed in Table 1 are relatively well determined, we can use them 
as a starting point to study these systems. 

A striking observation is that the eccentricity of at least one 
of the planets can be very high. Because hierarchical systems 
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90 



180 



270 



360 




180 270 
q) = 05 2 - wj (deg) 



360 



360 



Fig. 3. Possible secular trajectories for the HD 74156 system seen in the 
(o»i,ei) plane for two different values of the initial mutual inclination, 
Iq = 40° (top), and Iq = 50° (bottom). We show the trajectories us- 
ing the quadrupolar approximation, corresponding to the level curves 
of constant energy. Each one corresponds to an initial value of the ar- 
gument of the periastron aj t , ranging from 0° (black) to 90° (red) with 
a step of 10°. The dashed black curves corresponds to the observed ec- 
centricity of the planet, that is, it gives the initial condition for u>i . All 
trajectories are compatible with the current knowledge for this system 
(Table 1). 



exhibit high eccentricity values, the mutual inclinations can also 
be very high (e.g. Chatterjee et al. 2008). However, currently 
their mutual inclinations are unknown, not only because we are 
unable to determine the inclination with respect to the plane of 
the sky, /, (and hence the true masses), but also because we are 
unable to determine the longitude of the nodes in the plane of 
the sky, Q,: 



cos / = cos/i cos/2 + sin/i sin/2Cos(Q2 - Oi) 



(42) 



The only partial exception is the system HD 38529, for which 
h ~ 48°, estimated using astrometric measurements from the 
Hubble Space Telescope (Benedict et al. 2010). The planetary 
masses nt\ and m-i given in Table 1 correspond to the minimum 
masses (assuming I\ - h - 90°), except for HD 38529, where 
the masses were estimated using I\ - h - 48°. 

According to expression (3), without knowing the mutual in- 
clination of these systems it is impossible to determine the total 
angular momentum H, which is a constant of the motion. Even if 
we assume that the system is coplanar and prograde (cos / = 1), 
the total angular momentum is undetermined because the true 
planetary masses appearing in the expressions of G\ and G2 are 
also unknown (except for HD 38529). As a consequence, the 
present dynamics of these systems can be considerably differ- 
ent, depending on the true H value. 
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Fig. 4. Possible secular trajectories for the HD 74156 system seen in 
the (<p, e[) plane for / = 0° (top) and in the (o»i, e t ) plane for / = 35° 
(bottom). All trajectories are compatible with the current knowledge 
for this system (Table 1). We show the trajectories using the octupolar 
approximation (blue) and direct numerical simulations (red). In (a) the 
blue path corresponds to the level curves of constant energy for copla- 
nar orbits. The dot marks the present position of the planet. In (b) we 
additionally show the trajectories using the quadrupolar approximation 
(green), which corresponds to the level curves of constant energy. 



Moreover, although the argument of the periastron of these 
planets is known in the plane of the sky (angle <p* = - of 2 in 
Table 1), the arguments of the periastron in the invariable plane 
of the system (<x>\ and 0J2) are also unknown because they depend 
on /, and Q,- (e.g. Giuppone et al. 2012): 



cos(wi - of,) 



cos L - COS Ii COS / 

sin I\ sin / 



(43) 



where to* is the argument of the periastron in the plane of the sky. 
For a>2 we have an identical relation, the only difference is that 
a>i is measured from the ascending node, while u>2 is measured 
from the descending node. 

Therefore, apart from the semi-major axes, eccentricities and 
minimum masses, currently there are few constraints on the re- 
maining orbital parameters of hierarchical systems. Because we 
are only concerned with the tidal evolution of the mutual incli- 
nation, we adopted the masses listed in Table 1 for the numerical 
simulations. The only free parameters are then the mutual incli- 
nation and the arguments of the periastron, whose uncertainty 
are related to the lack of knowledge of the longitude of the node 
in the plane of the sky (O2 - ^i)- 

Using the quadrupolar approximation for gravitational inter- 
actions, the potential energy (Eq. 4) is independent of the a>2, and 
therefore e2 is constant (Eq. 10). The dynamics of the system is 
then fully described by the couples (/, u>i) or (e\,oj\). Adopting 
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the minimum masses (Table 1), we show in Figure 2(a,b) some 
possibilities for the HD 74156 system for different mutual incli- 
nations (corresponding to different H values). For nearly copla- 
nar systems (/ < 20°), only small variations are observed for e\ 
and /. However, as the inclination increases, the dynamics of the 
system is considerably perturbed by the presence of Lidov-Kozai 
cycles (Lidov 1961, 1962; Kozai 1962). In this regime, we can 
observe significant exchanges between the inclination and the 
eccentricity of the inner orbit. In some cases the eccentricity can 
reach values much higher than today, and thus enhance the tidal 
dissipation. 

Figure 2(a,b) is different from the standard Lidov-Kozai di- 
agrams that show level curves of the quadrupole Hamiltonian at 
fixed values G\ cos I, because for HD 74156 the initial eccen- 
tricity of the inner planet is already fixed at e\ = 0.64 (Table 1). 
Instead, we followed the procedure in Giuppone et al. (2012) and 
show the trajectories for different values of the initial mutual in- 
clination /, i.e., we varied the total angular momentum H. This 
explains why librating orbits in the Lidov-Kozai regime do not 
encircle a Lidov-Kozai equilibrium point (which occurs at the 
current eccentricity value). 

The impact of the initial uncertainty on the u)\ value is shown 
in Figure 3. Depending on this value, the observed eccentricity 
can correspond to a maximum or to a minimum for an identical 
total angular momentum H (Fig. 3a). Moreover, two trajectories 
may be in circulation or in libration (Fig. 3b). It is therefore very 
important to completely explore the phase space of the initial 
conditions (/, wi) to capture the global dynamics of hierarchical 
two-planet systems. 

3.3. Octupole contribution 

So far, we restricted our analysis to the quadrupolar gravitational 
interactions, because they are mainly responsible for the inclina- 
tion variations. However, for the hierarchical systems listed in 
Table 1, the range of semi-major ratios is 0.03 < a\la-i < 0.1, 
meaning that octupolar interactions cannot be neglected. Indeed, 
although the impact of octupolar terms on the eccentricity vari- 
ations is weaker than that of the quadrupole terms, octupolar 
interactions are strong enough to produce secular drifts when 
combined with tidal effects (Correia et al. 2012). 

In the planar prograde case (cos 1=1), the potential energy 
(Eq. 4) only depends on ip = a> 2 - oj\, and therefore / is constant 
(Eq. 12). The dynamics of the system is then fully described by 
the couples (<p, e\) or (<p, e 2 ). Adopting the minimum masses, the 
angle ip* in the plane of the sky listed in Table 1 corresponds to 
the angle (0)2 - o>\ + 180°) in the invariant plane of the system, 
because for I\ = I 2 = 90° we have u>\ = a>\ + 90° and w 2 = ^2 ~ 
90° (Eq. 43). In Figure 4(a) we show the expected eccentricity 
variations for the HD 74156 system in this unique situation for 
which the system is fully characterized. 

With increasing mutual inclination, we are left with four free 
parameters (e\,a>i, €2,0)2), and it becomes impossible to cap- 
ture the dynamics of the systems in a two-dimensional plot. 
Nevertheless, we can perform numerical simulations of the equa- 
tions of motion (Eqs. 13-16) to understand how the octupolar 
terms modify the quadrupolar approximation (Fig. 2). We ob- 
serve that the main effect is to add some diffusion around the 
quadrupolar trajectories. The diffusion is more pronounced for 
orbits in circulation around the separatrix (35° < I < 145°). The 
stability of the orbits can be measured with a frequency analysis 
(Laskar 1990, 1993). We determined the precession frequency g 
and g' of the argument of the pericenter o)\ over two consecutive 
time intervals of length T = 5 Myr. In Table 2 we compute the 



difference D = \g-g'\/g, which is a measure of the chaotic diffu- 
sion of the trajectories (Correia et al. 2005; Couetdic et al. 2010). 
It should be close to zero for a regular solution, and values with 
D > 10~ 6 correspond to chaotic motion. This is observed for all 
trajectories in circulation close to the separatrix. 

In Figure 4(b) we plot simultaneously the eccentricity evo- 
lution obtained with the two approximations for Io = 35°, and 
compare it with direct numerical simulations. We conclude that 
1) the quadrupolar approximation correctly describes the aver- 
age dynamics in inclined hierarchical systems; 2) the octupolar 
approximation is essential to derive a more realistic behavior and 
obtain results similar to direct numerical simulations. 



4. Tidal evolution 

We now include the effect of tides described in section 2.2 to 
the conservative equations of the motion (Eqs. 13-16), and per- 
form some numerical simulations. In all simulations we adopt 
for the innermost planet £ = 1/5, £2 = 1/2, and a dissipation 
time lag At = 200 s. For HD 74156 this dissipation is equivalent 
to k 2 /Q « 1.4 x 10~ 4 (Eq. 31), which is comparable to the value 
ki/Q = (2.3 ± 0.7) x 10~ 4 estimated for Saturn (Lainey et al. 
2012). In addition, we always set 102 = 180° and 2nl to = 50 day. 
The impact of the initial u> 2 value can be obtained by adjusting 
a different value for lo\. Similarly, the initial rotation rate is not 
a critical initial parameter, since tidal effects quickly bring the 
rotation near to the equilibrium value (Eq. 37). 

4.1. Effect of the spin 

In Figure 5 we show some examples for the evolution of the 
HD 74156 planetary system, in three different situations. The ra- 
dius of the inner planet is estimated to be R = 1 .23 R Jup (Eq. 41), 
and we initially assume / ( ) = 40° and a>\ = 0°. 

In a first experiment, we only consider tidal effects on the 
orbit (Eqs. 38, 39), as it is often done in previous studies. That 
is, since the rotation of close-in planets evolves very fast, we as- 
sume that the spin is locked in its equilibrium position (Eq. 37). 
We observe that the eccentricity and the semi-major axis of the 
inner orbit slowly decrease, while the mutual inclination and the 
eccentricity of the outer orbit only oscillate around a constant 
mean value (Fig. 5a). Therefore, we conclude that in this case the 
only effect of tides is to circularize the inner orbit in a timescale 
longer than the age of the system, a well-know result in the lit- 
erature (e.g. Correia & Laskar 2010b). 

In a second experiment, we neglect the effects on the orbit 
and we only consider the effect on the rotation (Eq. 25). This sit- 
uation corresponds to the opposite of the previous one, and it is 
not realistic, but it allows to highlight the importance of not ne- 
glecting the rotation rate evolution. Indeed, although there is no 
direct dissipative contribution to the eccentricity or to the incli- 
nation, we observe that these two parameters undergo significant 
variations (Fig. 5b), the eccentricity of the inner planet rising al- 
most up to 1 . The pumping effect on the eccentricity due to the 
spin excitation was reported in the planar case by Correia et al. 
(2012), for which only octupolar terms are important. In the non- 
planar case, the pumping effect it is even more pronounced, since 
quadrupole order terms additionally contribute. In addition, be- 
cause the angular momentum is mainly exchanged between the 
inner planet eccentricity and the inclination, while the first in- 
creases, the second decreases. In Appendix A we provide the 
full explanation for this effect in the frame of the quadrupolar 
approximation. We also observe that the eccentricity of the outer 
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Fig.5. Long-term evolution of the HD 74156 system for different tidal models with I = 40°: only dissipation in the orbit is considered (a), only 
dissipation in the spin is considered (b), full model (c). We plot the mutual inclination / (top), the eccentricities e\ (blue) and (green) [middle), 
and the semi-major axis a, (bottom). 



planet is slightly damped during this process, because of the oc- 
tupole order interactions (Correia et al. 2012). 

Finally, since orbital and spin evolution cannot be dissoci- 
ated, we integrate the full set of equations for the tidal evolution 
(Eqs. 25-28) (Fig. 5c). We observe that the initial behavior of 
the system is identical to the situation without dissipation on the 
orbit (Fig. 5b), with a significant damping of the mutual incli- 
nation. However, as the eccentricity increases, the inner planet 
comes closer to the star at periastron, and tidal effects on the or- 
bit become stronger. As a consequence, the semi-major axis de- 
creases and the damping effect on the eccentricity (Eq. 27) over- 
rides the pumping drift (Eq. A. 18). At this point, the inclination 
damping is less efficient, and it ceases when the pumping drift is 
completely gone. The system ultimately evolves into a circular 
orbit as usual, but in a considerable much shorter timescale and 
it is left with a final lower mutual inclination (~ 15°). 

4.2. Effect of the argument of periastron 

In Section 3.2 we have seen that the initial choice of the argu- 
ment of the periastron of the inner planet, o»i, plays an impor- 
tant role on how the present eccentricity is changing (Fig. 3). 
Similarly, it also changes the initial trend of the inclination: for 
increasing eccentricity e\ the inclination decreases, and vice- 



versa. In the example from previous section (Fig. 5c), we used 
oj\ = 0°, that is, we assumed that the observed value of the ec- 
centricity (e\ = 0.64) is a minimum, and the inclination / = 40° 
a maximum (Fig. 3a). 

In order to test the impact of the initial argument of the peri- 
astron in the HD 74156 system, in Figure 6 we plot its evolution 
for different initial to\ values. We observe similar behavior as 
before for all situations, the only significant difference being the 
evolution timescale. Until <x>\ < 50° this timescale increases, be- 
cause the present eccentricity is no longer a minimum value. As 
a consequence, the average value of the eccentricity oscillations 
is shifted down, and tidal dissipation is less effective, since at the 
periastron the inner planet is farther from the star. 

In the quadrupolar approximation, the behavior described 
above should be maintained up to lo\ < 90°, for which the ob- 
served eccentricity is a maximum and the inclination a minimum 
(Fig. 3a). However, the fact that the initial inclination increases 
when the eccentricity decreases has a strong implication when 
including octupolar terms (Fig. 2): for high inclinations the tra- 
jectories are closer to the separatrix, which results in a higher 
oscillation of the eccentricity. Thus, for o>i > 60° we observe 
that the evolution timescale is reduced again, since the inner or- 
bit eccentricity is allowed to reach much higher values than those 
predicted by the quadrupolar approximation. 
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Fig. 6. Long-term evolution of the HD 74156 system with I = 40° for different values of the argument of the periastron, a>i = 30° (a), 50° (b), and 
70° (c). We plot the mutual inclination / (top), the eccentricities e\ (blue) and e-i (green) (middle), and the semi-major axis a\ (bottom). 
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Fig. 8. Circularization time (e\ < 0.01) for the HD 74156 system with 
/o = 40° (red) and / = 50° (blue) using different values of the initial 
argument of the pericenter. For Io = 40° the timescale decreases around 
the Lidov-Kozai equilibria, because the system is the circulation regime. 
For Iq = 50°, the system is in libration, but it is most likely destroyed 
after it crosses the separatrix. 



Up to now, we have been considering an initial mutual in- 
clination /o = 40°. For Iq 55 40°, there is only one dynamical 



regime for the HD 74256 planetary system, consisting of trajec- 
tories in circulation around the Lidov-Kozai equilibira (Fig. 3a). 
However, for higher values of the initial inclination we can also 
observe the libration regime (Fig. 3b). In Figure 7 we show some 
numerical simulations with initial Iq = 50° using different values 
for the argument of the periastron. 

For u)\ = 30° the inner planet is still in circulation (Fig. 7a), 
so we observe identical behavior for the eccentricity and incli- 
nation as in the case with Iq = 40° (Fig. 6). At the bottom of 
Figure 7 we plot the inclination as a function of the argument of 
the periastron (/, u\). We plot a dot each 10 5 yr and its color be- 
comes darker with time. There we can clearly see that the planet 
is always in circulation, and that the amplitude of the oscillations 
is damped with time. 

In the remaining two situations shown in Figure 7, the inner 
planet is in libration around the Lidov-Kozai equilibrium located 
at <x>\ = 90°. The effect of tides is to slowly increase the ampli- 
tude of both the eccentricity and inclination. As a consequence, 
the orbit of the planet will cross the separatrix of the libration 
zone and start to circulate as in the previous case. 

For oji = 60° (Fig. 7b) the planet already starts close to the 
separatrix, so the initial oscillations are higher and it takes only 
about 75 Myr to cross it. For u)\ - 90° (Fig. 7c) the planet is 
placed close to the Lidov-Kozai equilibrium, so it takes much 
longer to reach the separatrix. In both situations, just after the 
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Fig. 7. Long-term evolution of the HD 74156 system with I = 50° for different values of the argument of the periastron, a>i = 30° (a), 60° (b), and 
90° (c). As a function of time, we plot the mutual inclination / (top) and the eccentricities e\ (blue) and ei (green) {middle). As a function of a>i, 
we plot again the mutual inclination, where the color of each dot becomes darker with time (bottom). 



transition of dynamical regime, the eccentricity reaches a very 
high value close to unity. Therefore, tidal effects with the cen- 
tral star become very strong and the final evolution is rapid: the 
semi-major axis decreases, and the inner orbit becomes circular. 
However, in a more realistic simulation where we integrate the 
full equations of motion and take into account the bodies dimen- 
sions, the planet most likely collides with the star. In both situ- 
ations the system is either destroyed, or its configuration com- 
pletely modified from the initial situation. 

In order to better understand the variation of the evolution 
timescale with the initial choice of the argument of the pericen- 
ter, in Figure 8 we plot the circularization time (e\ < 0.01) as 
a function of a>i . The circularization time is more or less equiv- 
alent to the inclination damping time, since the final stages of 
the evolution are very fast. For Iq = 40° the timescale decreases 
around the Lidov-Kozai equilibria, because the system is the cir- 
culation regime. For Iq = 50°, the system is in libration, but it is 
most likely destroyed after it crosses the separatrix. 

Giuppone et al. (2012) also studied the evolution of planets 
inside the circulation zone of Lidov-Kozai equilibriums. They 
performed some numerical simulations using the quadrupolar 
approximation and damping of the inner orbit eccentricity due to 
the presence of a primordial disk. They concluded that the planet 
stays in libration and migrates into the Lidov-Kozai equilibrium 



position, which is exactly the contrary that we observed here. 
Since the eccentricity of the inner orbit is also damped in our 
model, these results appear somehow contradictory. Therefore, 
we performed one simulation where the eccentricity is damped, 
but the semi-major axis is held constant. In this unrealistic situ- 
ation, one observe that the planet migrates into the equilibrium 
like in Giuppone et al. (2012). As a consequence, it seems that 
there is no inconsistency between the two models, but it becomes 
clear that the semi-major axis evolution plays an important role 
in destabilizing the Lidov-Kozai equilibria. It appears that it can- 
not be neglected in future studies on the migration of the initial 
orbits as in Giuppone et al. (2012). 

4.3. Constraints for the mutual inclination 

In previous sections, we saw that in mutually inclined hierar- 
chical two-planet systems there is a significant increase in the 
eccentricity of the inner planet's orbit. As a result, the tidal dis- 
sipation is enhanced when the planet is at the periastron, and the 
system evolves faster into an equilibrium configuration. 

In Figure 9 we show the evolution of the HD 74156 system 
for three different values of the initial inclination Iq = 15°, 30° 
and 45°. As expected, when we increase the mutual inclination, 
the evolution timescale decreases. For Iq = 15° the eccentric- 



10 



A.C.M. Correia et al.: Tidal damping of the mutual inclination 



ity and the semi-major take more than 10 Gyr to be completely 
damped, while for 7q =45° the system is fully evolved only af- 
ter 100 Myr. Moreover, for 7q = 15° there is almost no effect on 
the mutual inclination, we only observe some amplitude damp- 
ing when the eccentricity is decreased to low values, because 
gravitational perturbations no longer force the inclination. On 
the contrary, for 7o = 30°, the pumping effect on the eccentricity 
is already present, and hence we observe a significant reduction 
of the final mutual inclination. 

In Figure 10 we show the same kind of evolution as be- 
fore, but for an initial retrograde orbit with initial inclinations 
Io = 165°, 150° and 135°. In this case, the evolution of the sys- 
tem does not differ much from the prograde situation, the only 
significant difference is that the inclination is damped to high 
values close to 180° (coplanar system with a retrograde orbit). 
For /o = 165° the inclination is more or less conserved and the 
eccentricity is damped over 10 Gyr, while for lower values of 
the initial inclination, the inclination is damped and the system 
evolves in much shorter timescales. 
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Fig. 11. Circularization time (e t < 0.01) for the HD 74156 system us- 
ing different values of the initial mutual inclination and argument of 
the pericenter. Since the estimated age of the system is several Gyr 
(Table 1), we can rule out mutual inclinations within 20° < I < 150°. 



From the observation of Figures 9 and 10 we then conclude 
that mutual inclinations closer to 90° speed up the final evo- 
lution of the system. Since most hierarchical systems listed in 
Table 1 (except HD 190360 and HD 38529) still present substan- 
tial values for the inner-orbit's eccentricity, we then expect that 
their mutual inclinations are not extremely high. For all those 
systems we run several numerical simulations starting with the 
present initial conditions from Table 1, adopting k 2 At = 100 s, 
and different initial values for 7 and u>\. All trajectories that cir- 
cularize the inner-orbit in less than ~ 10 Gyr can then be ruled 
out, while those not showing significant modifications can be re- 
tained as possible representations of the real system. Therefore, 
we are able to set some constraints for the maximal mutual in- 
clination of each system, whose limitations are listed in Table 1. 
In Figure 11 we show the outcome of these simulations for 
the HD 74156 system, which corresponds to a summary of the 
more detailed evolutions shown in previous Figures. Orbits with 
20° < 7 < 150° circularize the system in less than 10 Gyr, so 
they can be discarded. 

When we run the same kind of simulations for the HD 38529, 
we observe that the eccentricity is damped very quickly, even 
for coplanar orbits. One possibility is that we overestimated the 
dissipation. However, even if we adopt k 2 At = 10 s, that is, one 
order of magnitude lower than for the remaining planets, the sys- 
tem still circularizes in a timescale shorter than the age of the 
system. Another possibility is to suppose that the inner planet 
semi-major axis was higher in the past. This hypothesis can also 
be extended to the HD 190360 system, for which the inner orbit 
is already circularized, but it may have had a higher eccentric- 
ity value in the past. In order to test this scenario, for all plan- 
etary systems in Table 1, we run several numerical simulations 
for different initial values for 7 and o)\, keeping k 2 At = 100 s, 
but adopting a\ = 0.2 AU and e\ = 0.25 as initial values, in- 
stead of the current values. In Figure 12 we plot an example 
for the HD 190360 system. By modifying the initial conditions, 
we are able to reproduce the present observations. Note that the 
HD 190360 system is older than the HD 38529 one (Table 1), 
so both systems may have undergone an identical evolution, but 
they are observed at different stages. The initial semi-major axis 
could have been higher, providing that the inner orbit eccentric- 
ity was also higher (for instance a\ - 0.25 AU and e\ = 0.4). 



4.4. Effect of the dissipation rate 

In the former sections we have been using k 2 At = 100 s, or, in 
terms of Q-factor, k 2 / Q ~ 1 .4 x 1 0~ 4 (Eq. 31), which is similar to 
the present value measured for Saturn kzl Q = (2.3 ± 0.7) x 10~ 4 
(Lainey et al. 2012). We adopted this value mainly for a better 
comparison with the previous paper on the planar case (Correia 
et al. 2012). Other works on the tidal evolution of hot-Jupiters 
also adopted similar values for At (e.g. Fabrycky & Tremaine 
2007; Correia etal. 2011). 

However, the Q-factor f planets is unknown, and may vary 
by some orders of magnitude (e.g. Goldreich & Soter 1966). 
Indeed, the value measured for Jupiter appears to differ by a fac- 
tor of ten, k 2 /Q = (1 .1 ± 0.2) x 10~ 5 (Lainey et al. 2009). On the 
other hand, statistical studies on the observed eccentricity distri- 
bution of hot-Jupiters give k 2 /Q ~ 10~ 6 (e.g. Jackson et al. 2008; 
Hansen 2010). Part of the problem is that the nature of tidal dis- 
sipation in these planets is still poorly understood. In addition, 
the Q-factor is frequency dependent (i.e., model dependent), and 
therefore sometimes it is difficult to translate from one system to 
another. In order to test the robustness of the inclination damp- 
ing, here we test our model with lower dissipation rates (higher 
Q values). 

According to the tidal equations (Eqs. 25-28), the evolu- 
tion timescales are linearly proportional to Q, so higher Q- 
values delay the final evolution of a system. In Figure 13 we 
plot the evolution of the HD 74156 system for three different 
dissipation values, k 2 At = 1, 10, 100 s (which is equivalent to 
kilQ ~ 10~ 6 , 10~ 5 , 10 -4 , respectively), starting with 7 = 50°. 
We observe that, although the evolution timescale is longer for 
higher Q values (as expected), the inclination damping is still 
present. 

In section 4.3 we saw that when we increase the initial mu- 
tual inclination 7o, the evolution timescale decreases very fast. 
For instance, in the case of the HD 74156 system, for 7o = 40° 
the inner orbit becomes fully damped after 1 Gyr, for 7 ( > = 50° it 
takes about 100 Myr, and for 7 = 60° only 10 Myr. Therefore, 
for different dissipation rates we can still put constraints on the 
mutual inclination of hierarchical systems, the only consequence 
is that as we increase Q, the maximal mutual inclination that one 
can expect to observe also increases. 
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5. Conclusion 

Many two-planet systems have been reported in hierarchical 
configurations. For most of these systems the mutual inclinations 
are unknown. However, since the orbital eccentricities are typi- 
cally high, we may expect that the formation mechanism that 
increased the eccentricities also increased the mutual inclina- 
tions. Very often the innermost planet in these systems is close 
enough to the star to undergo tidal dissipation, which can pose 
constraints on the final evolution. 

Here we have studied a particular subgroup of hierarchical 
planetary systems in which the inner planet's semi-major axis 
0.1 < fli < 0.5 AU (we called these "moderate close-in plan- 
ets"). This range is important to ensure that the spin of the inner 
planet is fully evolved, but not its orbit. Using an averaged sec- 
ular model that takes into account gravitational interactions up 
to octupole order, we showed that for many initial conditions the 
mutual inclination is damped to relatively low values (/ ~ 15° 
for HD 74156) on timescales shorter than the age of the system 
(less than one Gyr). 

Without planetary perturbations and for zero obliquity there 
is no effect from tides on the evolution of the inclination (Eq. 28). 
The inclination damping is thus not a direct consequence of tidal 
effects on the orbits. The key element is a inner planet's eccen- 
tricity oscillation at a secular timescale similar to the synchro- 
nization time of its spin. Indeed, the rotation (and thus the flat- 



tening (Eq. 2)) of the planet is driven by its eccentricity varia- 
tions (Eq. 25). In response to these excitations, the rotation is 
phase-shifted (Eq. A. 14) and the lag tends to pump the eccen- 
tricity (Eq. A. 18). Because the total angular momentum must be 
conserved, the increase in the inner orbit's eccentricity is ac- 
companied by a subsequent reduction of the mutual inclination. 
When the eccentricity pumping ceases, the inclination damping 
also stops. 

For high mutual inclination values, quadrupolar gravitational 
exchanges with the eccentricity are more efficient, and so is the 
eccentricity pumping. As a consequence, the inner orbit's eccen- 
tricity reaches higher values, tidal effects are enhanced at the pe- 
riastron, and the system evolves on shorter timescales. A strong 
inclination damping is then often followed by a fast circulariza- 
tion of the inner orbit. Since most of the observed hierarchical 
systems still present substantial values of the inner-orbit's eccen- 
tricity after several Gyr, we expect that they cannot have very 
high mutual inclinations. In particular, we are able to set con- 
straints on the highest mutual inclination in these systems. 

The evolution timescale also depends on the argument of the 
periastron in the invariant plane of the system, which is unknown 
at present for most systems. Indeed, the uncertainty on this pa- 
rameter is much higher than the uncertainty on the dissipation 
time lag Af. If a>i is in circulation, the evolution can be ex- 
tremely fast (a few million years) for values close to u)\ = 90° or 
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cj\ = 270°. However, for high mutual inclinations a>\ can be in 
Lidoz-Kozai libration for values close to cl>\ = 90° or cl>\ = 270°. 
In that case, the evolution timescale can be delayed to several 
Gyr, but the equilibrium is unstable and broken when the system 
crosses the separatrix with the circulation regime. After that, the 
inner planet is most likely lost. 



Appendix A: Inclination damping 

According to expression (28), the direct contribution of tidal ef- 
fects to the inclination can be neglected for low obliquities. We 
adopted e = 0° for simplicity, but in a more realistic situation 
the obliquity is expected to be trapped in a Cassini state (e.g. 
Correia et al. 201 1). Although high-obliquity states are possible, 
they are unlikely for close-in planets (Levrard et al. 2007), so we 
only expect these planets to be trapped in low-obliquity Cassini 
states, and hence dl/dt as 0. 

When tidal effects are combined with planetary perturba- 
tions, some counter-intuitive behavior may appear (e.g. Wu & 
Murray 2003; Mardling 2007; Correia et al. 2012). In particular, 
it has been shown that for moderate close-in planets, the gravita- 
tional perturbations of a distant coplanar companion combined 
with tides can increase the eccentricity of the inner orbit to very 
high values (Correia et al. 2012). Because of the total angular 



momentum conservation, the eccentricity of the outer orbit is si- 
multaneously decreased. 

For coplanar systems (sin / = 0), the octupolar terms in the 
expansion of the potential (Eq. 4) are the main responsible for the 
eccentricity variations of the inner orbit (Eq. 13). However, for 
mutually inclined systems, the quadrupolar terms become dom- 
inating, and we can neglect the octupole contribution: 

§(1 -e 2 ) 1/2 sin 2 / 

<?1 ~ V 2 l ^-rj- <?i SH^W] . (A.l) 

(1 -effi' 1 

In addition, the contributions to the outer orbit eccentricity van- 
ish (<?2 ~ 0), meaning that the angular momentum exchanges 
only occur between e\ and / (Eq. 3). Thus, if somehow the 
quadrupole interactions are able to pump the inner orbit eccen- 
tricity as in the planar case (Correia et al. 2012), one can expect 
to observe a decrease in the mutual inclination of the system. 

To understand this mechanism, we can simplify equations 
of motion without loss of generality. Since we assume that the 
orbits are not coplanar, we retain only the quadrupolar terms for 
the gravitational perturbations, that is, the conservative motion 
can be described solely by expression (A.l) and the first four 
terms in expression (15): 

VQ Vl x 2 

Ul * (l-e\) + (l-e^ 2 
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+ V21- 



+ v 22 - 



2(1 - 4) + | (ef - sin 2 /)(l - cos2wi) 
(1 - e 2 ) 1 ^ - e 2)3/2 

(1 + |e 2 - fe 2 cos2wi)cos/ 



2" 1 



(1 



e 2 ) 2 



(A.2) 



We can neglect tidal effects on orbital quantities 
(Eqs. 26, 27), which is justified since K' <k K (Eq.30). 
The only contribution of tides is then to the rotation rate 
(Eq.25). The semi-major axis and the mean motion are thus 
constant, as is the eccentricity of the outer orbit and G 2 . 

For simplicity, we average expression (A.2) over a)\, and lin- 
earize the set of equations of motion in the vicinity of the aver- 
aged values of x, e\, and I. Let x — xo + Sx, where xq is the 
solution of (37), e\ — eio + 5e\ , and I — Iq + 61. In the following, 





bO 




50 


\—>j 




CD 




T3, 


40 


C 




o 


30 


CO 








o 


20 


c 






10 









60 




oU 


o5 




CD 




2, 


40 


C 




o 






30 






inci 


20 




10 




U 




ou 




50 




CD 






40 






o 




"5 


30 






inci 


20 




10 








1 1 1 1 1 1 1 1 


k 2 


At= 100 s 















































-I 1 — I I I I 1 1 1 1 1 — I I I I 1 1 1 1 1 — I I I I 1 1 

k 2 At= 10 s 




_l I I I I I I I 



-I — I — I I I I I 



k P At= 1 s 




j i I i i ' i i 



10' 



10 10 J 

time (Gyr) 



Fig. 13. Possible past evolution of the HD 74156 system with / = 50°, 
using different dissipation rates. We plot the mutual inclination I for 
k 2 At = 100 s (top), k 2 At = 10 s (middle), and k 2 At = 1 s (bottom). 

61 is expressed as a function of 8e\ using the angular momentum 
conservation (Eq.3): 



61 



1 I G 



sin/o \G2 



+ COS Iq 



-'10 



1 -e 1 
1 e io 



6e\ 



(A.3) 



Then, the set of equations of motion (A.l, A.2, and 25) re- 
duces to 



6e\ = A sin2wi , 

= g + g x dx + g e Sei 
6x — —v x 6x + v e 6e\ , 
with 



A = v 2 



f(l-e? )" 2 sin 2 / 
(1 -e 2 F 2 



-ew 



(A.4) 
(A.5) 
(A.6) 

(A.7) 
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v x = £/i(e 10 ) , 

v e = -K(f[(e m )x - / 2 '(eio)) , 



(A. 10) 

(All) 
(A. 12) 



where /((e) = 15(e + 3e 3 /2 + e 5 /8)/(l - e 2 ) 11/2 , and / 2 '(e) = 
3(9e + 65e 3 /2 + 125e 5 /8 + 5e 7 /8)/(l - e 1 ) 1 . 

At first order, the precession of the periastron is constant 
co\ - g, and the eccentricity is simply given from expression 
(A.4)as 



5e\ = -Aecos(2gf + 2w 10 ) , 



(A.13) 



where Ae = A/2g, and oj\ — gt + <jd w . That is, the eccentricity 
e\ presents periodic variations around an equilibrium value eio, 
with amplitude Ae and frequency 2g. Since g x 6x,g e 5e\ <K g, the 
above solution for the eccentricity can be adopted as the zeroth- 
order solution of the system of equations (A.4-A.6). With this 
approximation, the equation of motion of 5x (A.6) becomes that 
of a driven harmonic oscillator whose steady state solution is 



6x — -Axcos(2gf + 2<jj\() — (f>) , 



(A. 14) 



or, combining the two products of periodic functions, 

g e Ae 

6ei = A sin(2gr + 2w 10 ) - — A — A sin(4gr + 4w 10 ) 

4g 

s ,-Ai g ,A 

_sx — A sm{4gt + 4wio -</>)+ ^Axsin^. (A. 17) 
4g 4g 

The two middle terms in the above equation can be ne- 
glected since they are periodic and have a very small amplitude 
(g x Ax, g e Ae <s g). However, the last term in sin (p is constant and 
it adds a small drift to the eccentricity, 



5e\ 



with Ax = v e Ae/ yjv 2 . + 4g 2 , and sin cf> = 2g/ -y/v 2 . + 4g 2 . The 
rotation rate thus presents an oscillation identical to the eccen- 
tricity (Eq.A.13), but with smaller amplitude and delayed by an 
angle <p (see Correia 201 1). Using the above expression in equa- 
tion (A. 5) and integrating, gives for the periastron 

g x Ax g«Ae 
coi = gt+tom — - — sm(2gt+2LO\o-(f>)- — — sin(2gf+wio) .(A. 15) 



2g 



2g 



Finally, substituting in expression (A.4) and using the approxi- 
mation g x Ax, g e Ae <§: g gives 

5e\ « A sin(2gf + 2oj\q) 



2g 
gxAx 
' 2g 



A sin(2gf + 2a)\o) cos(2gf + 2cjiq) 
A sin(2gf + 2u w - 0)cos(2gf + 2a> w ) , (A.16) 



gxA 
4g 



Ax sin (f> 



v e g x A 2 



4g(v 2 x + 4g 2 ) 



(A. 18) 



Note that the phase lag <p between the eccentricity (Eq. A.13) and 
the rotation variations (Eq. A. 14) is essential to obtain a drift 
on the eccentricity. From expressions (A.7) and (A.8) we have 
A ~ V21 ~ g, while from expressions (A. 11) and (A. 12) we get 
Vx ~ v e ~ K. Consequently, 



sin0 



and 



Ax = 



2g 



g 



VvJ+4? y/K 2 + 4g 2 ' 



v e A K 1 

, , ~ - cos (f> , 

2g Vv 2 + 4g 2 2 y/K 2 + 4g 2 2 



(A. 19) 



that is, 



5e\ 



16 



Vl 



sin 2(f> sin 2<p 



(A.20) 



(A.21) 



Thus, the drift vanishes when = or 90 degrees, i.e., for strong 
dissipation (K » g), where 5e\ ~ vxg/K, or for weak dissipation 
(K <s g), where 5e\ ~ v\K/g, respectively. The pumping effect 
on the eccentricity is then maximized when = 45°, for K ~ g, 
with 6e\ ~ vi. 

We therefore conclude that the quadrupolar perturbation of 
an inclined companion enhances the effect of the octupolar terms 
in the planar case (Correia et al. 2012). However, since the an- 
gular momentum is exchanged with the inclination, by assum- 
ing G2 » G\ we obtain from expression (A.3) that SI < for 
7 < 90°, and 61 > for 7 > 90°. That is, the mutual inclination 
is damped as long as the inner orbit eccentricity is pumped. 

The main difference when we consider the full non- 
linearized problem is that the drift in the eccentricity (Eq. A. 18) 
cannot grow indefinitely. Indeed, when the eccentricity reaches 
high values, the drift vanishes (Fig. 5b). Moreover, tidal effects 
are also enhanced for high eccentricities and counterbalance the 
drift (Eq. 27). Hence, the drift in the eccentricity is never perma- 
nent, neither is the inclination damping, although they can last 
for the entire age of the system. 
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